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We test the reweighting of the quark determinant of O(a) improved Wilson fermions in the 
domain-decomposed hybrid Monte-Carlo algorithm. Specifically, we implement a reweighting 
in a twisted-mass parameter proposed by Palombi and Luscher in Nf = 2 QCD. We find that 
at equal acceptance rate, the algorithm is significantly more stable on a 32 x 64 3 lattice upon 
switching on the reweighting parameter. At the same time, the reweighting factor does not fluctu- 
ate strongly and hence is under control. At equal statistics, the uncertainty on the pion correlator 
is comparable to the case of the standard, unreweighted algorithm. 
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1. Introduction 

Simulating Wilson type fermions at small quark masses is challenging, in particular because 
of the potential instabilities caused by the fluctuations of the low modes of the Dirac operator. A 
direct study of the low-lying spectrum of Wilson Dirac operators reveals that a spectral gap forms 
for large volume lattices [1]. However, occasional near-zero modes may appeal - since chiral sym- 
metry is explicitly broken by the discretization. In simulations based on the hybrid Monte-Carlo 
algorithm (HMC, [2]), this effect can lead to large 'spikes' in the history of the molecular dynamics 
Hamiltonian violation. Some time ago, Palombi and Luscher proposed [3] to reweight the fermion 
determinant so that the zero modes are suppressed by construction. The low-mode contribution 
to any observable is faithfully restored by including the reweighting factor in the ensemble aver- 
age. For a review of other applications of reweighting with various fermion discretizations, see the 
contribution of A. Hasenfratz at this conference. 

In the simplest version of the idea, a twisted mass term is added to the Dirac operator, 

D(ix)=D w + iiiY 5 , (1.1) 
where Dw is the 0{a) improved Wilson Dirac operator (see e.g. [4]). For Nf = 2, the weight 

W = det( D ™ Dw ) (1.2) 
\D w D w + p 2 J 

should be included in the expectation value of the observable G, 

(1.3) 



where (•••)„ stands for the ensemble average for the modified Dirac operator D(fi). 

Evaluating the reweighting factor W exactly is normally not possible, nor is it in fact required; 
instead it can be calculated stochastically. One may add a set of N pseudo-fermion fields (T]^, k = 
1, . . . , AO to the theory with action 

N 

^ = £(^,^), (1.4) 

k=l 

and the reweighting factor is replaced by 



W N 




D W D W + ^ 2 



dId 



w u w 




(1.5) 



The simulation with respect to the modified Dirac operator proceeds as before and the reweighting 
factor W is estimated, for each gauge configuration, according to Eq. (1.5) with N randomly chosen 
pseudofermion fields. Obviously, the larger the number of pseudofermion fields, the more accu- 
rate the estimate. Our tests show that 40 pseudofermion fields are sufficient to yield satisfactory 
reweighted measurements of meson correlators on 32 3 x 64 lattices. As is visible from Eq. (1.5), 
fluctuations in the reweighting factor W that strongly suppress the contribution of certain config- 
urations to the statistical average can occur if the original Dirac operator Dw admits very small 
eigenvalues. 
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Convincing arguments were given by Palombi and Liischer that the reweighting factor would 
not end up being exponentially small in the volume, as one would at first expect. This follows 
from the fact that the fluctuations of the lowest eigenvalues become smaller when the volume is 
increased [1] (possibly with the exception of a few). Nevertheless an explicit test is needed to 
ascertain that the twisted-mass reweighting idea works in practice, and this is the subject of this 
work. 



2. Implementation 

We use the domain decomposition preconditioning [5] of the hybrid Monte-Carlo algorithm 
(DDHMC). Our simulation code is based on Liischer' s open-source DDHMC code [6]. In this 
algorithm, the whole lattice is divided into check-board coloured blocks. Using Q. or £1* to denote 
the union of white or black blocks respectively, the Dirac operator assumes the form 

D (Da D dn \ 
\D da * D a * J 

where Da, Da* denotes the Dirac operator respectively on block £2, £1* with Dirichlet boundary 
condition, and Dga is the sum of all hopping terms from the exterior boundary dQ. of Q. to the 
boundary dQ.* of Q.*. The fermion determinant is factorized into local block parts and a global 
part. For Nf = 2 it reads 

detD^Dw = detR^R ■ Y[ detD A D A , (2.2) 

A 

where A runs through every block and 

R = 1 - P d a*D^D d aD^D d a* ■ (2.3) 



With even-odd preconditioning on every block determinant, Eq. (2.2) can be further written as 

det^-n (QlQee) fet (QlQoo) det (Q- 1 qY (q- 1 q) , (2.4) 

where 



Q = Qee-QeoQooQoe (2.5) 



and the site ordering has been chosen so that Da has the form 

Qee Q 
Qoe Q, 



YsDa = I n I • (2.6) 

loo 



The forces of the MD evolution steps are 

(co,F A ) = 2Re (Q-'Qee^Sco (^Qee) 0a) -2ReTr ( ^ + ^ ) , (2.7) 



Qee Qo 

(co,F R )=2Re(R- l ^ R ,8R-^ R ) , (2.8) 

where 0a is the pseudofermion fields supported on the even sites of block A and <p R is the pseud- 
ofermion field residing on the block boundaries. 
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The modifications required to introduce the twisted mass term are relatively benign. For the 
modified Dirac operator Eq. (1.1), the fermion determinant becomes 

detDin^D^) = det{D w - i^Ys)(D w + iixy 5 ) (2.9) 

where the reweighting parameter /I appears with opposite sign for up and down quarks. The forces 
take the form 

(w,F A ) = 2Re(^Q{^)-\Q ee + i^ A ,d a ^Q(^)-\Q ee + i^^ (2.10) 

™ „ fQeeSQee , QooSQ oo \ 

- 2ReTr — ^ + — = =■ , 

((0,F R ) = 2Re(R(ii)- l fa,S (0 R{ll)- l fa) , (2.11) 

where 

Qi.ll) = Q ee +ill - QeoiQoo + WY^Qoe , (2-12) 

R- l {n) = l-Pda*D{nr 1 D da *- (2-13) 

Note that to calculate the global force Fr, one needs to solve the full Dirac equation in every 
integration step. It is therefore important to solve the equation efficiently. The deflation acceler- 
ated [7, 8], Schwartz-alternating-procedure (SAP, [9]) preconditioned GCR algorithm is applied. 
Since the deflation subspace need not be exact, it is constructed by a relaxation process, starting 
from a set of random quark fields \jfh which are updated iteratively according to 

\l/,^"D w u ' ¥l , (2.14) 

until the condition ||D^i///|| < Ml/// is satisfied, where M is in the range of low eigenvalues of 
(D^Diy) 1 / 2 , and "D w l " is an iterative procedure that approximates the inverse of Dw- In the 
current implementation [6], it is done with the SAP. The operator Dw is used in the relaxation 
procedure, the deflation efficiency for the operator D{[i) is equally good, because the deflation 
subspace needs not to be exact and pL normally is small in practice. 

The little Dirac operator, i.e. restriction of the Dirac operator D(fi) to the deflation subspace, 
is then specified as matrix 

Au(fi) = (\ff k ,D(ji)yri) . (2.15) 

The Dirac equation D(n)y(x) = t](x) is separated into two equations by acting with projectors 
Pl(h) and 1 - Pl(h) from the left, where 

N 

P L {n)v(x) = £ D(ii) ¥k (x)A(ii) kl l ( ¥l , ¥ ). (2.16) 

3. Testing 

We have tested the reweighting method for two flavors of &{a) improved Wilson fermions on 
large lattices (24 3 x 48 and 32 3 x 64) with fine lattice spacings (respectively 0.07 and 0.08 fm). The 
masses of the lightest pseudoscalar mesons in these ensembles are approximately 360 and 300 MeV 
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configuration configuration 



Figure 1: Reweighting factor W on 24 3 x 48 lattices for ajX = 0.00569 (left) and 0.003 (right). 

respectively. Some details of the simulation, for instance the molecular dynamics trajectory length 
T, step size 8x, the number of configurations (ncfg) of the generated ensemble and acceptance rate 
of the HMC simulations are listed in table 1 . 

The reweighting parameter fi must be chosen with some care. Larger fi accelerates the algo- 
rithm but leads to larger fluctuations in the reweighting factor W. We have tested two values of afi, 
0.00569 and 0.003, on 24 3 x 48 lattices. The reweighting factors W are calculated according to Eq. 
(1.5). Their Monte-Carlo history is displayed in Fig. 1, after being renormalized so as to have an 
average value of one. 

For aji = 0.00569 (left plot), the reweighting factor W fluctuates strongly: about 20% config- 
urations receive a weight that is smaller than 0.01. By contrast, for a\i = 0.003 (right plot), the 
reweighting factor fluctuates moderately about the mean value. We may use the kurtosis 

K=^-3, (3.1) 
a 4 

where [in is the fourth central moment and a is the standard deviation, to quantify the distribution of 
the reweighting factor W. For a\i = 0.003, the kurtosis is K ~ —0.1 1, showing that the fluctuations 
of W is close to normal distribution; while for a\i = 0.00569, kurtosis rises to K ~ 12.8. Our choice 
of ji on the 32 3 x 64 lattice and the corresponding kurtosis is also listed in table 1 . 

In our simulations, we have tuned the integration step size of molecular trajectories to achieve 
an acceptance rate of « 80%. For comparison, we have simulated the same set of physical parame- 
ters with the standard method (jU = 0), using the same molecular trajectory length. We have tuned 
the acceptance rate to be similar, ~ 80%. In order to do so, the integration step size had to be 
roughly 10% smaller. Although both simulations thus ran at similar acceptance rates, we observe 
that the molecular dynamics trajectories are more stable when reweighting is used compared with 



Volume 


P 


a/fm 




ncfg 


T 


8z 


acc. rate 


a\l 


K 


24 3 x 48 


5.3 


0.07 


360 


140 


0.5 


0.028 


0.83 


0.003 


-0.11 


32 3 x 64 


5.2 


0.08 


300 


120 


2 


0.013 


0.78 


0.001 


-0.37 



Table 1: Parameters of the simulation, as explained in the text. The step size 8t is for updating the global 
force Fr. 
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trajectory trajectory 

Figure 2: History of Hamiltonian changes of molecular dynamics trajectories in HMC simulations on 32 3 x 
64 lattices. The reweighting method (left) is more stable than the standard method (right). 

the standard method. We plot the changes in Hamiltonian AH of molecular dynamics for each 
trajectory in Fig. 2 for simulations on the 32 3 x 64 lattice (reweighting method on the left and stan- 
dard method on the right). For the reweighting method, the history of AH is fairly stable, while for 
the standard method spikes show up more frequently; the highest spikes are up to three orders of 
magnitude higher compared to the reweighted algorithm. 

Using the reweighting method, we have calculated pion correlators and compared them with 
the pion correlators measured on the ensembles generated in the standard way. Practically one 
first measures the correlator in the standard way on the ensemble generated with the reweighting 
parameter /»; we refer to this correlator as the 'partially quenched' one. Then one 'corrects' it 
configuration by configuration with the reweighting factor (1.5). In Fig. 3, we plot the partially 
quenched pion correlator (gray) and the reweighted one (red) and compare them with the pion cor- 
relator obtained from the standard simulation (blue). They are compatible within the uncertainties, 
and have comparable statistical errors. 



4. Conclusion 

We have implemented a proposal by Palombi and Liischer to perform a simulation with a quark 
determinant that is protected from the fluctuations of the low-lying modes. We were able to find a 
reweighting parameter for which the stability of the DDHMC algorithm is significantly improved, 
and at the same time the reweighting factor is under control. As a test observable, we found that the 
correlator of the pseudoscalar density is well behaved under the reweighting (1.3). The behavior of 
other observables, such as the nucleon correlator, remains to be explored. 

In addition to studying the performance of the reweighted algorithm over more than 1000 
trajectories, it would be worth trying out other versions of reweighting, where the net benefit in 
stability is potentially even larger. Palombi and Liischer for instance made a proposal in this direc- 
tion [3] in which the UV modes are only affected at order /I 4 . 
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Figure 3: Pion correlator on the 32 3 x 64 lattice, with and without reweighting. 
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